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1. Introduction 

The origin of the large-scale structure is nowadays understood from the gravitational 
collapse of initial density perturbations which were produced by amplification of 
the quantum fluctuations in the inflaton field [T]. The properties of the large-scale 
structure depend both on the initial conditions at the end of inflation and on the 
growth of perturbations in a universe filled with non relativistic matter and radiation. 
The theory of cosmological perturbations is thus a cornerstone of our understanding 
of the large-scale structure. The evolution of radiation (photons and neutrinos) needs 
to be described by a Boltzmann equation [H |31 [Jj . Two types of perturbative schemes 
have extensively been used in the literature in order to describe the evolution of 
the cosmological perturbations. The first is a 1 -|- 3 covariant splitting of space-time 
m [7] and the second is a more pedestrian coordinate based approach. In the first 
approach, exact equations on the physical space-time are derived and perturbative 
solutions around a background solution are then calculated. In the second approach, 
an averaging procedure is implicitly assumed and, starting from a background space- 
time, perturbation variables satisfying the equations of motion order by order are 
constructed. In the 1-1-3 approach, the variables used are readily covariant, but 
the absence of background space-time can be a problem to simplify the resolution by 
performing a mode expansion, since the Helmholtz function is in general not defined on 
the physical space-time. In the coordinate based approach, all perturbation variables 
live on the background space-time, and enjoy the advantages of its highly symmetric 
properties. However, this extra mathematical structure is at the origin of the gauge 
issue through the identification mapping that we needs to be defined between the 
background space-time and the physical space-time. Thus, the gauge dependence 
needs to be understood. An elegant solution is to construct gauge-invariant variables 
a la Bardeen both for the metric perturbation variables ^ and for the distribution 
function j9l llOj. Since the Boltzmann and Einstein equations are covariant, they can 
be expressed solely in terms of gauge invariant variables provided we have a full set 
at hand. A full comparison of these two formalisms has been performed at first order 
in Ref. [11] , and for gravitational waves at second order in Ref. [12] . 

In the coordinate based approach, the true degrees of freedom identified from 
the Lagrangian formalism, are quantized. They transfer to classical perturbations 
which inherit a nearly scale invariant power spectrum and Gaussian statistics, when 
their wavelength stretches outside the horizon, thus providing initial conditions for the 
standard big-bang model. Conserved quantities [T3]fT4] enable to ignore the details of 
the transition between infiation and the standard big-bang model (see however |15j), 
and the evolution details need only to be known when the wavelength reenters the 
horizon. A first step to extend this procedure in the 1 + 3 formalism has been taken 
in Ref. [16] where conserved quantities were defined. As for the degrees of freedom 
which need to be quantized, a first proposal was made in Ref. [17], in order to identify 
them, but it has not yet been motivated by a Lagrangian formulation. 

The properties of the observed cosmic microwave background (CMB) anisotropics 
have confirmed the validity of the linear perturbation theory around a spatially 
homogeneous and isotropic universe and have set strong constraints on the origin of 
structures, as predicted by inflation. It now becomes necessary, with the forthcoming 
increasing precision of data that may allow to detect deviation from Gaussianity 
|18| . to study the second-order approximation, in order to discuss the accuracy of 
these first-order results. These non-Gaussian features are also of first importance. 
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since they can help discriminating between different inflation theories. Indeed, one- 
field driven inflation leads to very small levels of primordial non-Gaussianity [19j . 
whereas multifield inflation can present significant non-Gaussian features pn].|21j. 
However, since non-Gaussian effects also appear through non-linear evolution, that 
is from the second-order approximation and beyond of the evolution equations, the 
study of second-order evolution equations is necessary in order to distinguish between 
primordial and evolutionary non-Gaussianities (see Ref. |22j for a review on non- 
Gaussianity). Second-order Einstein and Boltzmann equations have been written in 
the 1-1-3 formalism [231 [21] i but not integrated numerically, partly because the mode 
expansion is not defined on the physical space-time, and this would then require a 
four dimensional numerical integration. However, the promising formalism of Ref. [25], 
which builds a bridge between the 1-1-3 formalism and the coordinate based approach, 
might shed some light on these issues. Similarly, in the coordinate based approach, 
the second-order Einstein equations have been written in terms of gauge-invariant 
variables [26] , and a first attempt has been made to write the Boltzmann equation in 
a given gauge for the different species filling the universe, and to solve them analytically 
[271 [28]. 

The goal of this paper is to provide the full mathematical framework for handling 
distribution functions at second order in the coordinate based approach taking into 
account the gauge issue. This will clarify the existing literature and point out some 
existing mistakes. We first review briefly in section II the gauge transformations and 
the procedure to build gauge invariant variables. We then present in section HI the 
transformation properties of the distribution function, and express them up to second 
order. We define in section IV the gauge-invariant distribution function and the gauge 
invariant brightness up to second order in the particular case of radiation (but this 
is readily extendable to cold dark matter). We then deduce in section V, from the 
Boltzmann equation, the evolution of the gauge invariant brightness in its simplest 
collisionless form, at first and second orders. To finish, we express in section VI the 
fluid limit as a consistency check of our results. 



2. Overview on gauge transformations and gauge- invariant variables 

2.1. First- and second- order perturbations 

We assume that, at lowest order, the universe is well described by a Friedmann- 
Lemaitre space-time (FL) with flat spatial sections. The most general form of the 
metric for an almost FL universe is 



ds2 = .g,,^da;^da;'' (1) 
= a{'qf{ -{1 + 2^)dTf + 2w,da;M?7 -I- [(1 - 2^')% -I- /i^^ ] dxMx^ } , 

where 77 is the conformal time and a the scale factor. We perform a scalar-vector-tensor 
decomposition as 

Lo, = d,B + B, , (2) 

/ly = 2Sy + d,E, + d,E, + 2d,djE, (3) 

where Bi, Ei and Eij are transverse {d'^Ei ~ d''Bi — d'^Eij = 0), and Eij is traceless 
= 0). There are four scalar degrees of freedom ($, ^E", B, E), four vector degrees of 
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freedom (S,;, Ei) and two tensor degrees of freedom {Eij). Each of these perturbation 
variables can be split in first and second-order parts as 

W = W^^^ + . (4) 

This expansion scheme wih refer, as we shah see, to the way gauge transformations 
and gauge-invariant (GI) variables are defined. First-order variables are solutions of 
first-order equations which have been extensively studied (see Ref. [29j for a review). 
Second-order equations will involve purely second-order terms, e.g. W^"^^ and terms 
quadratic in the first-order variables, e.g. [Ty*^^-*]^. There will thus never be any 
ambiguity about the order of perturbation variables involved as long as the order of 
the equation considered is known. Consequently, we will often omit to specify the 
order superscript when there is no risk of confusion. 

At first order, 4 of the 10 metric perturbations are gauge degrees of freedom and 
the 6 remaining degrees of freedom reduce to 2 scalars, 2 vectors and 2 tensors. The 
three types of perturbations decouple and can thus be treated separately. As long as 
no vector source terms are present, which is generally the case when no magnetic field 
or topological defect is taken into account, the vector modes decay as a~^. Thus, we 
can safely discard them and set i?^-^^ = B^^^ = 0. In the following of this work, we 
shall not include vector modes for the sake of clarity. We checked that our arguments 
and derivation can trivially (but at the expense of much lengthy expressions) take 
them into account. 

In the fluid description, we assume that the matter content of the universe can 
be described by a mixture of fluids. The four-velocity of each fluid is decomposed as 

u^ = i((5^+t;^). (5) 

The perturbation has only three independent degrees of freedom since must 
satisfy w^u^ — —1. The spatial components can be decomposed as 

v' = d'v + v' , (6) 

v"^ being the vector degree of freedom, and v the scalar degree of freedom. The stress- 
energy tensor of this fluid is of the form 

Tf_,i' = pu^Ui, + P {g^^ + Ufj,u^) , (7) 

where the density and pressure are expanded as follows 

p^p + Sp, P = P + dP. (8) 

At the background level, the form of the stress-energy tensor is completely fixed by 
the symmetry properties of the FL space-time. However, at the perturbation level, 
one must consider an anisotropic stress component, tt^^ with tt'^ = u'^vr^jy — 0. The 
pressure and density of the fluid are related by an equation of state, P — p/i, in the 
case of radiation. 

At first order, the formalism developed by the seminal work of Ref. [8] provides 
a full set of gauge-invariant variables (GIV). Thanks to the general covariance of the 
equations at hand (Einstein equations, conservation equations, Boltzmann equation), 
it was shown that it was possible to get first-order equations involving only these 
gauge-invariant variables. In addition, if these gauge invariant variables reduce, in a 
particular gauge, to the perturbation variables that we use in this particular gauge, 
then the computation of the equation can be simplified. Actually, we only need to 
compute the equations in this particular gauge, as long as it is completely fixed, and 
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then to promote by identification our perturbation variables to the gauge-invariant 
variables. Thus, provided we know this full set of gauge invariant variables, the 
apparent loss of generality by fixing the gauge in a calculation, is in fact just a way to 
simplify computations. Eventually we will reinterpret the equations as being satisfied 
by gauge invariant variables. The full set of first-order gauge-invariant variables is 
well known and is reviewed in Ref. [29] and Ref. [30]. As gauge transformations up 
to any order were developed, it remained uncertain |31j . whether or not a full set of 
gauge-invariant variables could be built for second and higher orders. This has been 
recently clarified [55], and the autosimilarity of the transformation rules for different 
orders can be used as a guide to build the gauge- invariant variables at any order. We 
present a summary of the ideas presented in Ref. [31] about gauge transformations 
and the construction of gauge- invariant variables [26j . 

2.2. Points identification on manifolds 

When working with perturbations, we consider two manifolds: a background manifold, 
A^Oj with associated metric g, which in our case is the FL space-time, and the physical 
space-time A^i with the metric g. Considering the variation of metric boils down to a 
comparison between tensor fields on distinct manifolds. Thus, in order to give a sense 
to "Sg{P) = g{P) — g{Py\ we need to identify the points P and P between these 
two manifolds and also to set up a procedure for comparing tensors. This will also be 
necessary for the comparison of any tensor field. 

One solution to this problem [31' is to consider an embedding 4 -f 1 dimensional 
manifold A/" = x [0, 1], endowed with the trivial differential structure induced, 
and the projections V\ on submanifolds with Vq{N) = x {0} A^o and 
Vi{N) ^ M X {1} = Ml. The collection of AIa = Vx{N) is a foliation of TV, 
and each element is diffeomorphic to the physical space-time AJ i and the background 
space-time A(o- The gauge choice on this stack of space-times is defined as a vector 
field X on M which satisfies X'^ = 1 (the component along the space-time slicing R). 
A vector field defines integral curves that are always tangent to the vector field itself, 
hence inducing a one parameter group of diffeomorphisms 0(A, .), also noted 0a(.), 
a flow, leading in our case from (i>{Q.p G Voi-^)) = p € Vo{Af) along the integral 
curves to e VoiAf)) = (? G Vi{N). Due to the never vanishing last component 
of AT, the integral curves will always be transverse to the stack of space-times and 
the points lying on the same integral curve, belonging to distinct space-times, will be 
identified. Additionally the property X'^ = 1 ensures that (I)\^x{'Pq{N)) — V\{N), i.e. 
the flow carries a space-time slice to another. This points identification is necessary 
when comparing tensors, but we already see that the arbitrariness in the choice of a 
gauge vector field X should not have physical meaning, and this is the well known 
gauge freedom. 

2.3. Tensors comparison and perturbations 

The induced transport, along the flow, of tensors living on the tangent bundle, is 
determined by the push-forward and the pull-back (f)\ [32] associated with an 
element <f>\ of the group of diffeomorphisms. These two functions encapsulate the 
transformation properties of the tangent and co-tangent spaces at each point and its 
image. Indeed, the pull-back can be linked to the local differential properties of the 
vector field embedded by the Lie derivatives along the vector field in a Taylor-like 
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$*^j(x'^)=:r^ + Ae + -e:tr + --. (10) 



fashion (see Ref. [32] or Ref. [21]) 

k — OC . 1^ 

'i>i,A(r) - E (9) 

k=0 

for any tensor T. 

A remark about coordinates changes is on order here. When the tensor T is a 
coordinate x'^ (once fj, is fixed, it is a scalar field), the previous definition reduces to 
the standard finite coordinates transformation. 

A2 
2 

This is the standard way of defining an active transformation on the manifold, 
by transporting a point of coordinates x'^ to a point of coordinates x''^. This 
transformation, when performed on the coordinate system - considering the 
coordinates as a grid on the manifold that one would displace according to the active 
transformation - induces a passive coordinates transformation, if we decide that the 
new coordinates of a point q are the coordinates of the point p such that (f>\{p) — q. 
When considering a transformation induced by a field ^, we will refer to the passive 
coordinates transformation induced by the active transportation of the coordinates 
system. 

The expansion of Eq. ([9]) on Va{Af) provides a way to compare a tensor field 
on V\{N) to the corresponding one on the background space-time Va{JV). The 
background value being To = £'yT|-pQ(jv), we obtain a natural definition for the tensor 
perturbation 

k=i 

The subscript X reminds the gauge dependence. We can read the n-th order 
perturbation as 



(12) 

which is consistent with the expansion of perturbation variables of the physical metric 
in Eq. ([J]), since the physical space-time is labeled by A = f. However, the fact 
that the intermediate space-time slices P\{Af) are labeled by A removes the absolute 
meaning of order by order perturbations, as it can be seen from Eq. (|f f p . The entire 
structure embedded by J\f is more than just a convenient construction and this will 
have important consequences in gauge changes as we will now detail. 

2.4- Gauge transformations and gauge invariance 

If we consider two gauge choices X and Y, a gauge transformation from X to y is 
defined as the diffcomorphism 

(l^x^Y^x = (0XA)"'(0ya), (13) 

and it induces a pull-back which carries the tensor Ax 7a, which is the perturbation 
in the gauge X, to AyTx, which is the perturbation in gauge Y since 
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= AyTx + To. (14) 

As demonstrated in Ref. [3T] this family (indexed by A) of gauge transformations fails 
to be a one parameter group due to the lack of the composition rule. It should be 
Taylor expanded using the so called knight-diffcormorphism along a sequence of vector 
fields ^i. For the three first orders, the expression of this knight-diffeomorphism is 



^Y,xm = rx^YAxiT) (15) 

+ ^iC^.,+SC^,C^,+Cl)<^*x,xiT). 

The vector fields ^i, ^2 and ^3 are related to the gauge vector fields X and Y 
by ^1 = y - X, 6 = [X,Y] and = [2X - Y, [X,Y]]. By substitution of the 
perturbation by its expression in Eq. (|lip . we identify order by order in A, and obtain 
the transformation rules for perturbations order by order. The first and second order 
transformation rules, on which we will focus our attention, are 



S^pT - 4')r = 2C^/^^T + {C^, + 4JTo. (16) 

The fact that we had to follow n integral curves of n distinct vector fields for n-th 
order perturbations is a characteristic of knight-diffeomorphisms. It arises from the 
fact that, for the whole differential structure of M to hold, gauge changes are a more 
general type of transformations than simple vector-field induced flows. Consequently, 
the Taylor-like expansion must be of a more general type. Indeed, for a given gauge 
change between X and Y, the family of gauge changes (j)x^Y.x labeled by A is not 
always a group in A, and this happens for instance if [X, Y] ^ (See Ref. [21] for a 
graphic intuition). Although we could, for a fixed A = Aq, find ^ such that Eq. p6)) 
takes a form like Eq. pT|) up to a given order, for instance by fixing Aq = 1, and 
choosing 

+ + ^ (^ + ^[6,6]), (17) 

this would mean that intermediate space-times are useless, and we would then ask 
Einstein equations to hold only for Vq {Af) and Vxo (A/") . This would lead to equations 
in the perturbation variables that mix different orders. The resulting solution, for 
second order and above, would be very difficult to find. 
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2.5. Perturbed Einstein equations 

Instead, we prefer to use this more complicated but cleaner knight-diffcomorphism 
(Eq. [ini) to change gauge. It keeps the differential structure built on TV and we 
additionally demand Einstein equations to be satisfied on each V\{JV). This can be 
used to differentiate Einstein equations to first order w.r.t A and take the limit A — > 
in order to get a set of equations that formally take the form £i[S'^^'> g, — 0. 

Once solved for the solutions of the first-order Einstein equation, we can differentiate 
twice the Einstein equation w.r.t A and get an equation of the type 

f2[<5^'\g,<5^'^r] =5[(5(i).g,(5(i)T], (18) 

where S stands for a source term quadratic in the first-order variables (see (I2j for a 
concrete example). 

We see that the decomposition of perturbation variables in the form given by 
Eq. ([4]) will trigger a similarity between the equations, i.e. £i and £2 have the same 
form. Purely second-order perturbation variables satisfy the same equation as first- 
order perturbation variables do, but with a source term. With known sources and 
known solutions to the homogenous equation, the Green function method enables us 
to solve, at least formally, the second-order equations, and by recursion at any order. 
To summarize, the Taylor expansion "taylorizes" the process for solving the equations 
by dividing tasks among orders, since Einstein equations are satisfied order by order. 

2.6. Gauge-invariant variables 

General covariance, i.e. the fact that physics should not depend on a particular choice 
of coordinates is an incentive to work with gauge-invariant quantities. As we notice 
from Eq. a tensor T is gauge- invariant up to n-th order if it satisfies C^S^pT = 
for any vector field ^ and any r < n, as can be deduced by recursion. A consequence 
of this strong condition is that a tensor is gauge-invariant up to order n if and only 
if To and all its perturbations of order lower than n either vanish, or are constant 
scalars, or are combinations of Kronecker deltas with constant coefficients. Einstein 
equation is of the form G — T = 0, and for this reason is totally gauge invariant. 
However, we cannot find non-trivial tensorial quantities (that is, different from G — T) 
gauge- invariant up to the order we intend to study perturbations, with which we could 
express the perturbed set of Einstein equations. 

Consequently, we will lower our goal and we will build, by combinations of 
perturbed tensorial quantities, gauge-invariant variables. These combinations will 
not be the perturbation of an underlying tensor. This method will prove to be very 
conclusive since a general procedure exists for perturbations around FL. Eventually 
we shall identify observables among these gauge-invariant variables and the fact that 
they are not the perturbation of a tensor will not matter. It has to be emphasized 
that the transformation rules of these combinations are not intrinsic and cannot 
be deduced directly from the knight-diffcomorphism since they are not tensorial 
quantities. Instead, we have to form the combination before and after the gauge 
change in order to deduce their transformation rules. 

We now summarize the standard way to build gauge-invariant variables. For 
simplicity we consider only the scalar part of the gauge transformations, since we will 
not consider vector modes in the metric and fluid perturbation variables (again, this 
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could be done, but would just obfuscate the explanations). In the following, we split 
as 

^o = rW, e; = 9'L(''), with r= 1,2. (19) 



2. 7. First-order gauge-invariant variables 

In the subsequent work we present the transformation rules of perturbed quantities in 
a simplified notation. Instead of writing Wy'' — w'^^ + / (^i, ..,Cr)j in order to state 
that the difference between the expression of the r-th order pertubed variable W in 
gauge Y and in gauge X is a function / of the knight-diffeomorphism fields ^i, ...,^r, 
we prefer to write ly^'"' — > + / (^i, .., S,r)- We remind that the expressions of the 
fields (e.) 

i<i<r necessary for the knight-diffeomorphism are expressed in function of 
the gauge fields X and Y [see below Eq. p^ ]. From the transformation rules we 
deduce that the first-order perturbations of the metric tensor ([1]) transform as 



$(1) ^ $(1) +T(i)' +WT(i) (20) 

(21) 

_ ^(1) _ T^T^i) (22) 
£;(!) ^ ijd) + L(i) (23) 

E^^E^^ (24) 
while the scalar quantities related to matter transform as 

5(1) P ^5(i)p + p'T(i) 

i;W ^t,(i)_L(i'' (25) 
^(i)7r*^ (26) 

where a prime denotes a derivative w.r.t conformal time 77, and where Ti. = a' / a. 

From now on, we shall refer to these first-order transformation rules defined by 
^1 as Tj, ($(!)), r5,(S(i'),... or simply r($(i'), r(B(i'), ... For instance T{¥^'^) = 

We first note that the first-order tensorial modes and the first-order anisotropic 
stress are automatically gauge invariant. For the other perturbation variables, which 
are not automatically gauge invariant, they are two ways to understand the procedure 
to build gauge-invariant combinations. The first point of view in building gauge- 
invariant variables consists in finding a way to get rid of the undesired transformation 
rule. To do so, we remark that the combinations B^^^ — i?^^^ and —E^^^ transform 
under a gauge change as B^^^ - E'^^y S'l' - E^'^y - T^^), -E^^^ -E'^^'> - L^^). 
We can use these combinations to add ad hoc compensating terms to and ^'f-'^' 
by defining 

$(1) = $(1) + _ pd)')' + n (pd) _ pd)') (27) 

,i,d)^vl/d)_^(pd)_£;(i)'y (28) 



and are now gauge invariant, by construction. This can also be understood, 
from a second point of view, as a gauge transformation for and '^^^^ towards 
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the Newtonian gauge (NG) [Tj defined by s!l!j^q decomposed in T^^^q = B^) — 
, l'^^j^q = —E''^\ which transforms the perturbation variables as 

B^^^ (29) 
e'-^^ (30) 

$(1) _ $(1) = ci,d)^ = $(1) + n - E^^y) + - i?(i)')'(31) 
vi/d) _ xpd) = M,d)^ ^ M/d) _ n (^d) _ i^d)') . (32) 

Similarly the gauge-invariant variables that would reduce to Sp, SP and v are 

<5(i)p = S^j^IjP = S^'^p + -p' - ' 
<5(i)p = ^i^jP = ^(DP + P' (p(i) - £(1) 

7r'^(i) = 4^^7r'J' = <5(iV*^ (33) 

Since we have ignored the vector gauge degrees of freedom, i?^^ and i?*^^^ are 
the two gauge variant variables of the metric perturbation while and are 
the gauge-invariant part. As mentionned before, we then force the gauge-invariant 
variables in the perturbed metric by replacing <I>*^^) with <l'^^) — H. ^pd) _ ^ -|_ 

(ijd)_i5(i)') 

and applying similar procedures for ^I^d)^ jd)^ and (5*^^ -'P. When 

developping Einstein equations, we know that general covariance will eventually keep 
only gauge-invariant terms. Thus, we can either do a full calculation and witness 
the terms involving the degrees of freedom P'-^-' and P'-^-' disappear, or perform the 
calculations with pd) and P^' set to zero and obtain the perturbed Einstein equations 
only in function of gauge-invariant variables. The latter simplifies the computation, 
which is useful when going to higher orders. The advantage of the second point of view, 
is that the addition of the compensating terms of the first point of view can be seen 
as a first-order gauge change towards the Newtonian gauge with ^^^g. (decomposed 

as pI^tvg ^"^^ ^^-^ng)- These enables us to decompose the perturbed metric in a 
gauge-invariant part and a gauge variant part as 

<5(i)g = ,5(i).g + £ a, g, (34) 
as it can be seen from the transformation rules under a gauge change characterised 

by?! 

-^NG^ -^NG+^l- (35) 

This property which is not general but happens to hold in the case of cosmological 
perturbation (i.e. around FL metric) is the key to extend this construction to second 
order. 

It should be noted that this procedure, although achieved by defining gauge 
invariant variables which reduce to the perturbation variables in the Newtonian 
gauge, can be extended to other types of gauge-invariant variables which reduce to 
perturbation variables in another gauge. For instance, we can use the transformation 
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properties of and E'^^^ to add the compensating terms to S^^' and other 

variables. The transformation rules ^^^V^ ^ 'l!^^^'^ /H. -T'^^\ -E^^'> -E^^'^ - L'^) 
make it straightforward to build these compensating terms. We need to define C^^q 
decomposed with T^p^j = '^^^^ /H, L^^p^ = —E^^\ The gauge-invariant variables 
defined with this procedure reduce to the perturbation variables in the flat gauge 
= 0, ^-(1) = 0) and are 



B^Bfg = B^'^-^-EW^ $(i) ^ $W = ci>(i)+v|/(i)+f .(36) 
2.8. Second- order gauge-invariant variables 

For second-order perturbations, Eq. (jl6p gives the following transformation rules 

$(2) _ $(2) _,_ J^'(2) + ^y(2) ^ 

^(2) ^ ^(2) _ y(2) ^ ^,(2) ^ 

^(2) ^ ^(2) _ -^j.(2) _^ 

^(2) ^ ij(2) + i(2) + 

Eg^ -.E^f+SE., 

d^^^P ^ S'^^^P + P'T^^^ + Sp 

- 2tt'*'^^^ dkd' l''^^ - 2TT^^''^^dkd'L^^\ (37) 

where the source terms are quadratic in the first-order variables T^^^ ^'(^^ 
We collect the expressions of these terms in [Appendix A[ In the rest of this paper, we 
shall refer to these second-order transformation rules associated with (i^) = (^1,^2) as 
T(^) ($(2) ) 7-^^^ (5(2) ) ^ ^ ^ oj. simply r ($(2) ) ^ 7-(5(2) ) ^ ^ ^ ^ ^ These transformation rules are 

much more complicated than their first-order counterparts. However, the combination 
defined by F = 5^"^^ g + 2Cai) (5'^'<? + ^^i, g enjoys the simple transformation rule 

F F + Cg. , r,(i) , ,g under a gauge change defined by ^2 and ^1 (see Ref. 126 ). 
As a result, its transformation rule mimics the one of first-order pertubations under 
a gauge change. This means that if we decompose F in the same way as we did for 
the metric the metric with 

^<i>(2)+5^(^W^J 

Bp ^ b('^ + Spi^^c) 
Ep ^E(^^+Sp{^%a) 

Ep,, ^ E^f + SE,,{^%a), (38) 

then the transformation rules for these quantities will be similar to those of Eq. ((29|) , 
but with the vector ^2 + [C^wGi Ci] instead of ^1. Consequently, we shall use the same 
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combinations in order to construct gauge-invariant variables which are 

$(2) = $^ + {Bp - E'p)' + H{Bf- E'p) 

= '^F -TiiBp - E'p) 

E^f^Ep,,. (39) 

As for the first order, this addition of compensating terms can be understood, 
from the second point of view, as a defining the gauge-invariant variables as the 
perturbation variables in a given gauge. In our case it is the Newtonian gauge since it 

(2) 

transforms B and E into a null value. This transformation is defined by ^\}^q that 
we decompose in 



rp(2) 



L%a= - E^'' - S^^^ {^%a) ■ (40) 
The second-order gauge-invariant variables can be seen as 

E^^ EE S^^^E,, (41) 

r^^^'^ = s'-^ln'^ . (42) 



where the index NG indicates that we transformed the quantity with the formula ([T 
with the vector fields C^^^vg ^^jvg defined above. This means that we have split 
the second-order metric according to 

5(2)5^^(2)3 + /: g + 2£ a) 6^^^g^£\,r, g (43) 

where S'-^^^g is the gauge-invariant part and — ^^'jy^ the gauge variant part, as it can 
be seen from the transformation rules under a gauge change characterised by (^1,^2) 

5(2)5 -<5(2)~^ 

-^%G^ -e^L+6 + K^!lvG>Ci]. (44) 
As for the first order, we can choose other types of combinations, for instance 

those which are equivalent to setting the gauge as being flat, by using this procedure. 

(2) 

In this case, the vector field ^_^pQ is decomposed in 

T^a - ^^^S^' (^g) , L^o - -E^'^-S'^ (^'^o) -(45) 

It should also be mentioned that the existence of an inverse Laplacian A^^ of the 
background space-time, i.e. a corresponding Green function with boundary conditions, 
is required for all this procedure. In other words, when working in Fourier space, all 
our conclusions will be valid only for modes which do not belong to the Kernel of A. 
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3. Gauge transformation of the distribution function 

3.1. pre-Riemannian distribution Junction 

So far, we have set up the mathematical framework to identify points between the 
background space-time and the perturbed space-times through a gauge field X. This 
enabled us to define the perturbation of tensors and to calculate their transformation 
properties under a gauge transformation. However this only allows to perform a 
fluid treatment of the radiation. In the statistical description for a set of particles, 
we assume that each particle has a given impulsion and is located at a given 
position [33]. The equations then have to describe the phase space distribution of 
the particles. If the number of particles is high enough, we can define a probability 
density, the distribution function, of finding a particle in an infinitesimal volume of 
the phase space. Now, let us focus our attention on this distribution function. The 
distribution function is a function of the point considered (i.e. its coordinates a;^), 
and also a function of the tangent space at this point whose coordinate we label by 
p'^d^. There is no special reason for this function to be linear in p'^d^, but we can 
expand it, without any loss of generality, in power series of tensors according to 

(46) 

k 

The distribution function is then decomposed as the sum of all the multipoles J^^^,,^^ 
evaluated in a particular point of the tangent space. From the previous section we know 
the transformation rules for these tensorial quantities, thus / transforms according to 

r(j) [/(x^p'^)] ^ (47) 

where 7(j) refers to the knight-diffeomorphism with the set of vectors (^i, ^2, •■•)• 

As we do not necessarily want to refer explicitly to the decomposition in 
multipoles, we use the fact that for any vector ^ — ^''9^^, which defines a flow on the 
background space-time Vo{J\f), we can define an induced flow (a natural lift) on the 

vector tangent bundle TVoiM) directed by the vector field = ^^df,,?" {d,.^'')-^ • 
This implies the useful property 

(^mi..mJp^^-P^^ = ^T? {T^,..^^p^^..p''-) . (48) 
With this definition, we can rewrite the transformation rule for / as 

7^5)[/(x^p^)] = 7^T^)[/(^^p'')]: (49) 

where now TItq refers to the knight-diffeomorphism with the set of vectors 
(Ta,Te2,...)- ^ 

The evolution of the distribution function is dictated by the Boltzmann equation 
= C[/], where the r.h.s is the collision term which encodes the local physics, 
is collision term can be easily expressed in the local Minkowskian frame defined 
by a tetrad fields Ca, from known particles physics. For this reason, the framework 
developed to define gauge transformations for a general manifold has to be extended to 
the case of Riemannian manifold. Instead of using the coordinates basis 5^ to express 
a vector of tangent space as V = p'^d^, we prefer to use the tetrads basis ea and 
write V = 7r°ea. In terms of coordinates, this means that the distribution function 
is a function of and tt'^. When expressing the physics with the tetrad fields, the 
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metric is not just one of the many tensors of the theory whose properties under a 
gauge transformation we need to know, but rather a central feature of the manifold, 
since it determinates the tetrads (up to a Lorentz tranformation) required to express 
the distribution function. As the metric is a tensor, and as the tetrads are defined 
according to the metric, the extension is inherited from the previous section. 



3.2. Tetrads 

3.2.1. Definitions On each slice V\{J\f), we should have four vector fieldfJH (and 
their associated 1-forni fields) labeled by a = 0, 1, 2, 3, which satisfy the normalization 
conditions 

e^^e^<7p. = ^afc, e^ets'^'' = (50) 

With these notations, indices a, b, c. are raised and lowered with rjab- 

With the formalism developed for tensors, we carry this tetrad field onto the 
background space-time using a gauge field X with 



fc— oo . ^ 



k=Q 



An) ^ rn 

1 — ^X^a 



VoiJ^y e-. = 4°^e., (51) 



and similar formulas for e". 

As Ca is a basis of the tangent space on the background space-time (and e° a 
basis of its dual space), and eJJ can be expressed in the generic form 



where. 



e-a,X — Ra,X^b, ^X — ^""^a.Xi Ra,X^c,X — ^a,xRc,X — ^ ai (^2) 



Rab,X =^TiR', 



fc! "'''^ 

k 

fc! 



Sab^x^J2ljSal]x- (53) 



k 

Order by order, this reads 

dx ea--ti^,x^b, e - e b^ ^- (54) 

3.2.2. Normalization condition Tetrads are four vector fields which satisfy Eq. (|50p 
and are thus related to the metric. Consequently, the perturbations of the tetrad 
defined above are partly related to the perturbations of the metric. When pulled back 
to the background space-time, Eq. ([501) implies 



4'Xx{Vab) = Vab = (l>Xx«49t^:^) 

= ^lxK)4>lx{^'i)4>lx{9,u). (55) 

X The fifth direction which arises from the extension of the manifold from M to M is ignored as the 
component of any tensor is required to vanish in this direction. We thus consider the tangent space 
at each point of A/" as being four-dimensional. 
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Identifying order by order (in terms of A) we get in particular for the first and second 
orders 

eb.(5^''eQ + ea.(5^^eb + (5^''g(ea,eb) =0 

eb-S'^x^ea + ea.^^$eh + (5^''g(ea, e;,) (56) 



■'(5^^eb.(5^''ea + (5^-'.g [s''^ ea, + 5^^^ g {ea, (5^''e6^ 



= 0, 



n(2) 
^(ab),X 



where a dot product stands for 5 (_, _). From the constraints (j56p . we can determine 
the symmetric part of r!^^^ as 

<l).x= - \5x9{ea,e,) (57) 

= - \&f9{ea, -eb) - 4^5 (^2^6-, e-,) 

-4'^ff(e-a,<xe-^)-4'].^<x, (58) 
which are related to the components of the inverse by 

(59) 

■^ab.X - - KkX + ^K,xKb,X- (OL)) 

The antisymmetric part, i?[ab],xj still remains to be chosen as it corresponds to 
the Lorentz transformation freedom (boost and rotation), which is allowed by the 
definition ([50]) . A first and easy choice would be ^["^j x ~ ^ ^^^^ However, as 
mentioned above, we eventually want to decompose a vector p^d^ on tangent space 
as 

= ^"^a = TT'^e^d^, (61) 

and identify tt" with the energy and tt' with the momentum (although conserved 
quantities are generally ill-defined in general relativity, energy and momentum can be 
defined when performing perturbations around a maximally symmetric background 
[55] as it is the case here). When working with coordinates, we want to express 
physical quantities, as measured by comoving observers, i.e. observers of constant 
spatial coordinates, whose motion is defined by the 1-form (d?])^ [IB]. We thus 
require (e*^)^ ^ {drfj^j^, which is equivalent to choose >S'^"q ^ = for any n, where 
tti = 1,2,3. This choice allows us to fix the boost in S'^"-' by imposing the condition 
^[a'!o],x = ~'^[o"al],x = ''^talo),x- ^q. ^ implies that for any n 

{p + q = n, 
P> 1, 9> 1} 

it can be checked by recursion that this implies 



',,],X 

?(») 



We also fix the rotation by requiring >5'f".''^ .i = 0, and it can be checked similarly, by 



recursion on Eq. (p^ . that this implies R]^^.^ x~^- 
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3. 3. Gauge transformation of tetrads 

Under a gauge transformation, we can deduce the transformation properties of the 
tetrad from those of the perturbed metric. For simpUcity, we restricted to scalar and 
tensor perturbations, but this is completely general and can be easily extended to 
include vectors. In the FL case, we use a natural background tetrad associated to 
Cartesian coordinates eo — (9,,) /a, eb^ — [di) /a, in order to evaluate Eq. (l57|) . The 
notation hi refers to Lorentz (S0(l,3)) indices running from 1 to 3, whereas i is a 
coordinate index running from 1 to 3. When uselessly obfuscating the explanation, we 
will not make the distinction and change hi for i. We report the detailed expressions 
for the transformation of the tetrads for the first and second orders in [Appendix B[ 

4. Distribution function 

Now that the transformation properties of the tetrads are known, we turn to the 
general transformation of a distribution function /(a;'',7r°). 

^■l. Multipolar expansion 

Any function /(a;^,7r") can be expanded in symmetric trace free multipoles as [37j 

= (64) 



with 



(65) 



We do not need any additional identification procedure for the tangent spaces 
through a gauge field, in order to identify points of the tangent space of the slices 
TVx{N). Indeed, once the metric and a gauge field X are chosen, there exists a natural 
identification with the tetrad fields. First, and as mentioned before, we identify the 
points of A/" which lie on the same integral curves of X, that is, we identify a point 
P G Voi-N') and $a,x(-P) G V\{N). Then, we identify vectors of their respective 
tangent spaces, if the coordinates of these vectors in their respective local tetrad 
frames e-a and Ca, are the same. To be short, we identify Tr^Ca and 7r"ea. As a 
consequence, for any given set {ai,...,ap}, the function Fa-t^,,ap{x'^) is a scalar field. 
Fai..aj,{x'^) is then pulled back on the background space-time using the gauge field X, 
and we define in this way perturbations 



A" 



(66) 



and 



(67) 



This perturbation scheme induces a perturbation procedure for the distribution 
function / as 
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n 

p 

It is essential to stress that tt" is not a perturbed quantity, it is a coordinate of the 
locally Minkowskian tangent space. However, the tetrad field allows us to see as 
a perturbed vector since p''(7r°) = e^7r°. In other words, for a given tt", there is an 
associated vector whose order by order perturbation in a given gauge X is given by 

— „M(n)„a 



4-2. Gauge transformation: general case 



We can deduce the transformation rule under a gauge change directly on the form (|65|) , 
pulled back to the background space-time, 

p 

The first factor in this expression is tensorial. Exactly as for the pre-Riemannian 
case, its transformation rule is dictated by the knight-diffeomorphism, whereas we 
get the transformation rules of the tetrads from Eqs. (|B.2p and Eqs. (jB.4p . As we 
do not necessarily want to refer explicitly to the multipole expansion, the first factor 
is rewritten by considering / as a function of using 7r° = e'^xP'^^ ^^'^ ^^PPlyiJ^g 
Eq. ((49|) . Wc then have to consider the resulting distribution function as a function of 
7r°, knowing that the inversion is now given by p'^{tt°-) = T(e^)7r°. This will account 

for r f e 



in Eq. (j69p . In a compact form it reads 
r[/x(:.^^'^)]=r(To{/x[x^eX]} 



(70) 



To obtain an order by order formula, we explicit these three steps using a Taylor 
expansion. First, we use that 

d 



/x(a;^7^'^) = 



exp e->^VxTr^ /- 



IX 



(x^e>^)s<?x(x^p''),(71) 



in order to consider / as a function of p'^. We then Taylor expand back the result of 
the knight-diffeomorphism in order to read the result as a function of tt" , 



r[/x(x^^'^) 



d 



exp(et7r'^r(i?A)^)r(Te) (gx) 



(x^e->'^).(72) 



The derivatives in the previous expressions have to be ordered on the right in each 
term of the expansion in power series of the exponential. When identifying order by 
order, we need to take into account the expansion in Rab and Sab, in the exponentials 
and also in the knight-diffeomorphism. 

We have provided the general transformation rules for the distribution function 
and we will specify now the transformation properties of the first- and second-order 
distribution function. 
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4-3. The mass shell 

The transformation properties of (5^''e(^ have been chosen so that, in the special 
case of / = gf^uP'^'p" — g^„e^^e'j^TT°-n^ — nair"', it remains unchanged under a gauge 
transformation, i.e. T{Tr'^TTa) = Tr'^Tr^. Since the tetrads must satisfy Eq. ([SD]) . then 
^x^f ~ ^ ^^or ri > 1, and it imphes this property triviaUy. As a consequence, any 
function of the form S{TraTr'^ — rn^)/(x^, tt") transforms as 5{TTaTT°- — m^)T [/(x^, tt"')], 
where is the mass of the particles described by the distribution function. In other 
words, the transformation of the distribution function remains on the mass shell, as 
it has been already mentionned in Ref. [9]. We will make use of this property when 
computing the transformation rules of the distribution function. 



5. Application to the perturbation of the Boltzmann equation for 
radiation 

The formalism developed in the previous section is general. We will now apply it to 
the particular FL case, and from now on we will also focus on the radiation case, 
that is the case where m? = 0. For the first and the second order, we will present 
the transformation rules of the distribution function for radiation, and build a gauge- 
invariant distribution function as well as a gauge-invariant brightness. We will then 
write the evolution equation of this gauge-invariant brightness in the case where the 
photon travels freely through space-time without being affected by diffusion processes. 
This is obtained using the coUisionless Boltzmann equation 

^ = ^ + ^^ + ^^ + ^^^0 (73) 
dT] drj dx^ drj dn'^ drj dn^ drj 

where n* = ttYtt", from which we will extract the background, the first- and the 
second-order equations after having pulled it back to the background space-time. In 
order to do so, we need to know and By considering p'^ as a perturbed vector, 
as mentionned in § (|4.ip . theses can be expressed from the geodesic equation 

P°^ - -K.P^P" (74) 

that we pull back to the background space-time in order to extract order by order 
equations. Similarly, is given by the order by order expressions of p"^- = P% 
when pulled back to the background space-time. 

At the background level, space is homogeneous and isotropic. Consequently, 
the distribution function depends neither on the direction of the photon nor on 
the position in space . It only depends on tt" and r/, which implies that = 

^ = 0. Since the background geodesic deviation equation implies — — TYtt", the 
coUisionless Boltzmann equation reads at the background level 



dl 
drj 



-H.°|^ = 0. (75) 



5.1. Gauge transformation at first order 

In order to better understand the seemingly heavy but powerful formalism of § 14.21 let 
us apply it to the first-order gauge transformation of the photon distribution function 
/ in the Boltzmann equation. In this case, Eq. ([70]) for = (T, L) leads to 
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r 



The expressions of Rj^x ^^"^ '^'^j^^^, and their transformation rules for the FL case, 
are given in [Appendix B[ Using the fact that / is only a function of 7r° due to the 
term 5{-KcT^'^)^ 



B - Bf 
[Kx^^ap^^)] = T— fix'', ap^^) + T^7r"(T' + n'B.T) 



T R, 



>(i)o 



(1)0 



BttO 



Btt^ 



T R 



Note that there is no term involving 
prescription in the choice of the tetrad in § 13.2.21 
We then express the derivatives as 



(1)0 
i.X 



0(1)0 
'^i,X 



af{x\ap^^) 



Bt] 



Btt° 



Putting all the pieces together, we finally get that 



(77) 



(78) 



^^'§^ thanks to the 



(79) 



T 



5(^e^^)4'7 ='5(7r%c) 



= (5(7r%c) ( -^TT^n'B^T + 



{Btt^ 



Br] 



(80) 



= <5(7r%,)^7r"(HT + n'9,r), 

where in the last step we have made use of the background Boltzmann equation (j75p . 

It can be checked that by considering / as a function of instead of 7r°, as 

allowed by the factor 5{'k'^'Kc), we recover the same result as performed in Ref. [9]. 
However this is slightly more intricate, as it now apparently depends on the three 
variables tt* which are in fact not independent at the background level. 

Although the mathematical framework can seem to be heavy, we did not need 
to define an extension of the distribution function outside the mass shell nor a gauge 
transformation field parallel to the mass shell as in Ref. We first have built the 
distribution function using the tetrad field (it is a function of 7r° and not an express 
function of p^). Then, as explained in ^14. 31 the normalization condition ||5D|) . when 
expressed at each order in Eqs. (j56p . ensures that it remains on the mass shell during 
a gauge transformation that we perform using the rules derived for tensors. 



5.2. First-order gauge-invariant distribution Junction for radiation 

Now that transformation properties of the first-order distribution function are known, 
we can use the results of § [2] to define a gauge-invariant distribution function by 



/(I) 



^(1) f 



Of 



n 



_ E^^y) + n'B, - E^^y)] . (81) 
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As for tensorial quantities, we can choose for instance S^L^pQ in the above 
expression, in order to define an other gauge-invariant distribution function. Its 
expression is given by 



5/ n 



n 



These two first-order gauge-invariant distribution functions are related by 



/(I) _ JUJ = 



(1) _ 



5/, 



4,(1) 



(82) 



(83) 



It is worth remarking that in the previous hterature , another gauge- invariant 
distribution is defined, namely 



^ f(i) + ^Tr^^-^i) 



(84) 



Though it cannot be interpreted as the perturbation of the distribution function in a 
given gauge since it mixes £,^ng ^-nd C^fGj this is a better variable to highlight the 
conformal invariance of the photons propagation and to compare with the null cone 
integration method 

This first-order analysis illustrates the power of this formalism which can be 
generalized to higher orders in perturbations. 



5.3. First-order collisionless Boltzmann equation for radiation 

Integrating the gauge- invariant distribution function of radiation over tt'', we define 
the gauge-invariant brightness, which is the energy perturbation per unit solid angle 
in a given direction 

i(i)(a;^,nO = 47r J f'^^l {xf" ,Tr'> ,n'){7T°fd7T° . (85) 

We choose the normalization of the background distribution function such that 
the background brightness reduces to the energy density (see § [5] for the fluid 
approximation) 



1(7^) = At: / /(,7,7r-)(vr"rd^" = p. 



(86) 



We can associate gauge-invariant symmetric trace-free moments, to this 

brightness by using the decomposition 



iW(x^n»)^^^^),^( 



x'^ ..n ■ 



(87) 



With these definitions, the integral / (tt^)'^ d7r" on the first-order Boltzmann equation 
leads to the evolution equation for I^^^ [35] 

(1^ + ^ + + - ¥'y) 1=0, 
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where we have ignored the tensor terms for simphcity. Similarly, a gauge-invariant 
brightness X^^-' associated with f'^^\ and a gauge-invariant brightness M'^^^ [S] 
associated with T'^^^ can be defined. They are related to l'^^^ by 

M(i' =±(1) -4X4'(i). (89) 
5.^. Gauge transformation at second order 

At second order, the general gauge transformation of the distribution function ([70|) 
for (0 = (6,^2), [TO = (r^i,T^2) is given in details in |Appendix~Cl After 
simplifications, it reads 



r ((5^V) = §^(r^^^ + TT' + d^Td'L) 



+ ^7r"|n*a,T(^) - 2n^ [(a.^ji; -|- + d.djL) d'T ~ *9jT] 
-f- d{Td'T + {Tn'd,T)' + n'^i {d^ LdjT) + 2^n'diT^ 

{7:°f{n'd,Tn^d,T)+2j^Tn'd,T+^T^ 



+ 2^hJ^n"n^djT + 2^§^7rOa^r + 2d'Ld/^^ f + 2r^^. (90) 

This is a cornerstone expression in our study of the second-order distribution 
function. As for the fluid quantities, knowing the transformation rules under a second- 
order gauge change is enough to define a second-order gauge invariant distribution 
function which is required to write the second-order Boltzmann equation only in 
terms of gauge-invariant variables. As for tensors, several gauge-invariant distribution 
function can be defined, and this relation is also required to express how the different 
gauge-invariant distribution functions are related. 

5.5. Second- order gauge-invariant distribution function for radiation 



Again, we can use the results of ? l2.8l to define a gauge-invariant distribution function 
as 

^ 4%/ = sff + r. ^ Uff] . (91) 

As for tensorial quantities, we can choose for instance (^C^^fg^ ^^to) ' ™ order 
to build another second-order gauge-invariant distribution function. 



?(2) _ r(2) r _ r(2) ^ , ^ 

\^ — * FG '' ^ FG I 



5ff . (92) 



The difference between these two gauge-invariant distribution functions is also 
gauge-invariant and is consequently expressed only in terms of gauge invariant 
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quantities. For the sake of completeness, we give the form of the relation between 
these two gauge-invariant distribution functions, 



zn \ 



2n V 



2 d'f 



n'd. 



J_dl 



■^^(2) + (1)' + 2n¥^'^'^ + n¥^ 



'#(1 



¥^^n'^^¥^') 



n 



2H 



( 9' 9,4(1)9^ a,*(i) 
2n \ ^ 

This clearly demonstrates the power of our formalism since, contrary to the first 
order, this relation cannot be guessed intuitively. Note also that this is non-local as 
it is generally the case for second-order gauge-invariant quantities. 

5.6. The second-order gauge-invariant collisionless Boltzmann equation for radiation 
Similarly to the first order case, we define the second-order brightness as 

i'^'^\x^' ,n') = All f /(2)(x^,7r°,nO(7r")3d7r". (94) 



We also define the second-order gauge-invariant moments associated to this gauge 
invariant brightness by the second-order version of Eq. ([87]) • The derivation of the 
collisionless Boltzmann equation in the Newtonian gauge is detailed in Ref. [571 HH] ■ 
Once the integral / (7r°)'^d7r*' performed, it leads to an evolution equation for the 
brightness. As this is a scalar equation, it is gauge invariant and it can be expressed 
only in terms of the gauge invariant quantities that we have defined and which reduce 
to the perturbation variables in the Newtonian gauge. Explicitly, it reads 



2X(1) (^(1)' -nJ9j$(i)) -l(xp(2)' +4^p(i)*«' 
l(¥^^ +¥'An'd,I^'^ =0. 



(95) 



Up to this stage, we agree with the expressions of Ref. [27l [28] . 
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6. Fluid approximation 

If we want to recover the transformation rule and the gauge-invariant variables for 
the energy density, the pressure and the velocity of radiation, we need to define a 
stress-energy tensor from the distribution function of radiation. We already know 
from special relativity how to define such a tensor. We generalize it by using the local 
Minkowskian frame 



T°^(a;'^) = / d7r°d37r',5(7r%e)/(a;^,7r^)7r"7r^ 

{-K°ff{x^',^T'^)n''n''A^:°d^n\ (96) 



where n°' = = ttYtt^, if a = 1, 2, 3 and n° = 1 if a = 0. In order to evaluate the 
stress energy tensor, we have performed one of the integrals which removes the Dirac 
contribution 5{-K°--Ka) 

(5(7r''7r„)G(a;,7r")d7r°d37r* = /G(x,7r°,n*)7r°d7r"d2n\ (97) 



Several useful relations for handling integrals of the background distribution function 
are reported in [Appendix D[ If we are dealing with several species, we can still define 
a stress-energy tensor for each component, as long as we are dealing with weakly 
interacting gases. This is the standard kinetic approach in which the interaction 
between two species is encoded in the collision term of the Boltzmann equation [2] . 
We define p, P, the velocity t/" and the anisotropic stress 11"'', 

T'^ = pV'U^ + P Jl"^ +^''^ (98) 

with l.°-^= 77"'' + WU'', and the properties VUa = -1, W'' ±ab= 0, UaU"'' = 0. 
However, fluid quantities are usually expressed using the canonical basis associated 
with coordinates and not the tetrad field. We thus define -u'' = U°-e'^ as 
the coordinates of the velocity in this canonical basis, and we decompose it as in 
Eq. ([5]). Similarly, we define the anisotropic stress expressed in the canonical basis 
by tt'^" = ej^ej^n"''. Some confusion can arise from the fact that physicists often 
design a vector by its coordinates. With this symbolic convention, U"" and are 
mathematically the same vector, but expressed in different basis since If^Ca = u^^df^. 
The relations between [/" and up to second order are 

C7" = as" = 1 

W = au' = 0, (99) 



and 



=5' fw(i) , (100) 



5fu^^d^{v + B)d'{v + B) 

(2) 



8fu' = a'(w(2) + b(2)) _ 2$a*B + 2^d' {B - v) 

+ 2d^v~B){d'djE + E'j). (101) 
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Similarly the relations between the spatial components of tt^" and XI"'' are 



, ^-(102) 

The fluid quantities can be extracted from Eq. as follows 

p - T'^'f/aC/fc, (103) 
3p ^ yafc ^^^^ (^04) 

Hafc = T^'^ (^±ca^dfc ^cd^ab^ , (105) 

{p + P)U°U' = T°\ (106) 

It is easy to see that the factor S (TTaTr") in the integral of the definition (|96p of the 
stress energy tensor implies that P — p/3. 

The system of definitions (|103m06| determines the fluid quantities. Indeed, these 
quantities can now be calculated itcratively at any order once Eq. ()96p is pulled back 
to the background space-time. Since U'^ = 1 and = 0, p P and 11°'' are given by 

p^3P = T™C7of7o, n"'' = 0, (107) 

as expected from the background symmetries for a fluid of radiation. Then, since 
U'^ = y/V^Ui + 1, and using Eq. (|106p . we can determine the first-order expression of 
the velocity 

Repeating this procedure, we obtain from Eqs. ()103m06p 
-34^)p = <5«r<'°c7oC7o 

= - —S^x^Tl , (109) 



and the condition = implies 

^'^noo = 

SPu°' =26''x'^U'^S^x^Uj. (110) 
Again, using Eq. (jl06p . we determine the second-order perturbation of the velocity 

5(^){7° = 4'^C/*4'^C/, (111) 
42)^. ^ ^ 1^42)^0, _ ^(2)^0,-) _ 2^4)L7^ (112) 
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Iterating, we obtain from Eqs (|103m06)) 

=34'^P = 4'^rO"c7oC/o + 2TO"c7o5|^C/o (113) 

This shows that, by iterating this procedure, the fluid quantities can be 
determined up to order n if /, that is T""^, is known up to order n. This means 
that, by knowing the transformation rule of / under a gauge transformation, we can 
deduce the transformation rules of the fluid quantities built out of it (p, P, [/°, 11'^''). 
Eventually, we are interested in their expressions in the canonical basis in order to 
compare with the results of § [H and we need to use Eqs. (|99m0ip and Eqs. p02p . 



6.1. First- order fluid quantities transformation 

At first order, from the relations (|109p and (|100p . and the transformation rule for /, 
Eq. ([50]) . we deduce after some algebra, that s'^x^ p transforms as in Eq. ([^5]) . Similarly, 
from Eq. (|108p . the relation (jlOOp . and the transformation rule for /, Eq. ([SO]) , we 
deduce that w^^^ transforms as in Eq. ([25]). By the same method, we recover easily 
that 6'-^'^ TT*-' is gauge invariant. 



6.2. First- order fluid equations 

In order to recover the gauge-invariant conservation equation and the Euler equation 
of the fluid approximation at first order, we define the first-order gauge invariant 
stress-energy tensor by 



4tt 

and its associated first-order gauge-invariant fiuid quantities, p'^^'> , P'-'^^ , w'^^^ and n'^^^^') , 
built from the same types of relation as in the set of Eqs. (| 10311 106)) and expressed in 
the canonical basis with Eqs (|100p and (|102p . Because of the comparison performed 
in the previous section, these quantities match those defined in Eq. (j33p . and this 
justifies the fact that we use the same notation. We need the useful relations between 
the first moments and the fluid quantities 

■pH) ^ fjiD^^S^'^p, (115) 
J 47r 

[iWn^^^-pd^v^^\ (116) 
J 47r 3 

P^^'^ = J i(i) (^nV - ^ = (117) 

Performing J dil on the brightness evolution equation (|88p . we recover the flrst-order 
conservation equation. However, performing J nMfi, we recover the first order Euler 
equation as expressed in Appendix E only if we neglect the first-order anisotropic 



pressure. This comes from the fact that the statistical description of radiation leads 
to an infinite hierarchy of equations coupling moments of order p — 1, p and p -I- 1 [39j , 
whereas the fluid description keeps only the equations involving the monopole and the 
dipole. 
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6.3. Second- order fluid quantities transformation 

In order to establish the second-order comparison with the fluid description, we need 
to know how to perform an integral involving S)^ f, for instance on 2 ir'^n^djT + 
2 Q^f TT^d'^T. We will thus make use of the multipolar expansion 

/ ^ f^P. + 4/a, + 5(1)) + ^4i)n,,nV^" + ... (118) 

from which it can be checked that we recover the correct fluid quantities when used 
to compute J^^^T"*' in Eq. (|96| . 

Using the same method as for the first order, with the relations (|113p and (|10ip . 
and the transformation rule for the second-order distribution function, Eq. ([90|) . we 
deduce that S^^^^ p transforms as in Eq. (P7|). Additionally, from the relations (lllip . 
(fm]) . pm]) and pIO)) . we deduce that v^^^ transforms as in Eq. ((37)l . 

We also notice that from the definition (|113p . the relations (|llll) (|112p . and the 
transformation rule for /, Eq. (PO)) . we deduce that ^'^^II*-' transforms according to 

j(2)n'J ^ 5(2)ny + 2T (^(i^ff^y -I- 2d''Ldk (^(i^ff^) . (119) 

When expressed in the canonical basis (tt^" = e^^e^II"''), we recover exactly the 
transformation rule of the anisotropic stress given in Eq. ([37|) . 

This is one of the major results of this paper. We recover the perfect fiuid 
transformation rules for the energy density, the pressure, the velocity and the 
anisotropic stress given in Ref. [40j up to second order, when starting from the 
statistical description. 



6.4. Second- order fluid equations 

In order to recover the gauge-invariant conservation equation and the Euler equation 
of the fiuid approximation at the second order, we follow the same procedure as for 
the first order case. We thus define the second-order gauge invariant stress-energy 
tensor by 

f"^(2)(^P) = y(^O)3/(')(x^7^)n'^n''d7r"d2^7-yi(2)„V^,(120) 

and its associated second-order gauge-invariant fiuid quantities, p^^^ P*^^) u'^^ and 
7j-»j(2)^ built from the same types of relations as in the set of Eqs. (|103m06|) and 
expressed in the canonical basis with Eqs pOip and ()102p . Because of the comparison 
performed in the previous section, these quantities match those defined in Eq. (f42|) . 
thus justifying the fact that we use the same notation. 

In order to recover the conservation and Euler equations of the fiuid 
approximation we perform the integral J ^ and / on this equation. However, at 
the second order this has to be done with care since the link between the second-order 
gauge-invariant brightness and the second-order fluid quantities is given by 

^(2) ^ f j(2)^ ^ + (121) 

J in 3 



p{2) ^ 



/ X(2)„'^ ^ -pf9^{;(2) _2#(i)a*{}(i)U-5(i)p9^z)(i),(122) 
J 47r 3 V / 3 
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\ dVL 
T j 4^ 



3 



(123) 



This clearly differs from the expressions (5.10) and (6.33) of Ref. [77] where the term 
quadratic in v in .F*-^-* , the term quadratic in ^' and v in JF*(^) are not there. The 
difference in the energy density perturbation as extracted from JF^^^ , comes from the 
fact that the fractional energy density A*^^^ for the radiation is defined as seen by the 
observer of velocity ~ (d??)^ whereas we define it in the fluid frame. The fractional 

energy density that they define is related to our quantities by pA^^' = (S^^T™f/oC/o. 
The difference in the expressions for the fractional energy density can be traced 
using Eqs. pi3p with Ea. ([llip . However, this is only a matter of definition and it 
is consistent with Eq.(7.2) of Ref. [27]. Implicitly the authors of Ref. [57] do also 
use a tetrad basis in their section 3 in order to identify coordinates of the tangent 
space between the background and the perturbed space-time, in the same way as 
explained below Eq. (j65p . Their p is equal to our tt*^ and the unit vectors match 
when restricting to the Newtonian gauge. The equations (3.6) and (3.7) of Ref. [27] are 
equivalent to Eq. (|6ip when expressed in the newtonian gauge with the use of Eqs. (|52p . 
Eqs. (|B.ip and Eqs. (|B.3p . As for the difference in the velocity perturbation as defined 
from ^ 

comes from the fact that their definition for w^'^' has to be interpreted in 
the tetrad basis, and therefore it matches Sj^qW. However, the difference between the 
tetrad basis and the canonical basis is not computed as in Eq. (|112p . and it explains 
the discrepancy. This can also be checked on the second-order extraction of Eq.(7.3) 
in Ref. [27] . Indeed, there is an the extra term quadratic in VE" and when compared 
to Eq.(2.15) of Ref. [H], as a trace of the difference between our perturbed velocity, 
which matches the definition in the canonical basis usually given by Eq.® and Eq.Q, 
and their perturbed velocity. However, the equations involving v'j^'^ in Refs. [571 155] 
such as Eq.(4.6) are consistent with this difference, though the physical interpretation 

i(2) 

Vj as being the perturbed velocity of photons in the canonical basis is not correct. 

The computation of a term like^^, is easily performed using the multipolar 
expansion 

fii) ^ _^ 4M«(i)n'^ -t- i^ni^^nV^" + ... (124) 

P 2p J 

Applying this method, we recover the second-order conservation equation detailed 
in [Appendix E As for the Euler equation, we recover it at second order only if we 



neglect the anisotropic stress up to second order (beware that the anisotropic stress is 
different from the second moment of the distribution as it can be seen on Eq. (|123p ). 
and use the first-order Euler equation. 

This is also a major result of this paper. We recover the fluid gauge invariant 
equations up to second order, only if we can neglect the anisotropic stress up to second 
order. It remains to be shown that this is extended up to any order, as we expect. 

Let us also stress that in Ref. ^8\, the term -§-r is evaluated using -§-r = 

^1^, in order to derive Eq.(4.1) and Eq.(4.6). However, this is not correct since / 
is a function of the independent variables r;, a;*, tt", n'. Even though they are related on 
a photon geodesic, they are independent in the analytic expression of /. Additionally 



Gauge-invariant Boltzmann equation and the fluid limit 



28 



this method is not fruitful because ^1 iv ~ Vi): since does not parameterize a 

photon geodesic. Consequently, the subsequent analytic expressions of this reference 
solving the conservation and Euler equation are not correct (for both radiation and cold 
dark matter) though the Boltzmann equation is correct. This can also be seen directly 
from the fact that these equations do not match fluid approximation equations of 
[Appendix E[ Once corrected for this mistake, and taking into account the differences 
mentioned before we can check that the collisionless part of the conservation and Euler 
equations (4.1) and (4.6) of Ref. [57| match our equations. 



6. 5. Validity of the fluid approximation in the literature 

In this paper, we have considered so far the fluid approximation as a theoretical 
framework in which we restrict the description of a species to its energy density and 
its velocity. The computations involved for the distribution function at second order 
were rather long, and it was used as a consistency check for the gauge transformation 
rules and the collisionless Boltzmann equation. Since the fluid approximation is built 
out of the kinetic theory, it was indeed expected that all the conclusions made in this 
statistical description could find their fluid approximation counterpart. 

It is now necessary to determine under which conditions this can be done, that 
is when the anisotropic stress can be neglected. This requires to work on the physics 
of coupled species, baryons and photons, in the cosmological context. The collision 
term as well as its physical implications have been studied in Ref. |28j and it is very 
likely that the extraction of its quadrupolc in Eq.(4.18) is not affected by the previous 
considerations. Indeed, in the tight coupling limit (which requires only the collision 
term) for a system of photons and electrons highly coupled through the Compton 
diffusion, the authors of Ref. [55] find that the quadrupole satisfies 



(125) 



This result is necessary to determine in which case the fluid approximation can be used. 
Comparing it with Eq. (jl23p . we immediately see that the physical interpretation of 
this result is that the second-order anisotropic stress of radiation W^^^^ is suppressed 
in the tight coupling limit. As a consequence, the fluid approximation can be used in 
the tight coupling limit also at second order in perturbations. 



7. Conclusion 



In this article, we have performed a general investigation of the gauge invariance of the 
distribution function. This allows us to recover very easily the standard results at the 
first order and to extend them at the second order. We derived the fluid approximation 
at first and second orders. This required to carefully define the stress-energy tensor 
in the local Minkowskian frame. At the second order, our results differ from the 
ones previously derived in the literature [27l |28] . We have tackled down the origin of 
the differences and shown that it was lying in an incorrect identification between the 
tetrad and the canonical basis. Our analysis, restricted to the collisionless case, puts 
the second order Boltzmann equation, needed if we intend to study non-Gaussianities 
in the CMB, on firm ground. 



Gauge-invariant Boltzmann equation and the fluid limit 



29 



Acknowledgements 

I thank Jean-Philippe Uzan for drawing the topic to my attention and for his endless 
comments on the manuscript. The second-order expansions were computed using the 
tensor calculus package xAct [4? and I thank Guillaume Faye and Jose Martin Garcia 
for their help on it. I thank Ruth Durrer Nicola Bartolo Sabino Matarrese and Antonio 
Riotto for commenting on their works. Finally, I thank Thiago dos Santos Pereira for 
his numerous remarks on the draft. 



Appendix A. Sources terms in second order transformations 



The perturbation variables in the decomposition ([T]) are extracted as follows 

1 

2a 
1 



* = - 77Z2^9Q0, (A.l) 



B ^^A-^^'5go^, 



Using this method we can read the source terms defined in Eq. (|37|). which 
are quadratic in the gauge change variables T, L and the perturbation variables 
^,^,B,E,E,j, in Eq. ^ 



S^^T (T" + 57iT' + [n' + 2n^)T + 47i$ + 2$') 
+ T' (2T' + 4$) + d,Ld' (T' + nT + 2$) 
+ d^L'd' {T-2B- L') , 

s^= - T {UT' + in' + 2n'^)T - 2^' - m.'H) 

- di [UT - 2*) &'L 
1 



(A.2) 



dj {2B + L' -T) diT 



did^L {2dkdjL + Adud^E + iEkj + {2nT ~ m)bkj) 
Td.dj (L' + 2nL) 

T {2E[. + 2d,djE' + AHEij + AHd^djE) 



d'^Ldk {d,d,L + 2E,j + 2d^djE) 



(A.3) 



is slightly different from Ref. [H] and Ref. [T3] since, in these works, the extraction 
of metric perturbation variables is not performed according to Eq. (jA.ip . However, 
this mistake does not matter for their study that focused on the long wavelength limit. 



Sb ^ A~^d'[T'di{2B + L' -T) 
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+ d^L' [2didjL + 2 (HT ~ 2^-) 5^ + 4 {E,j + d.djE)] 
+ d^d^Ldj (25 + L' -T) + Ldjdr {2B + L' -T) 
+ a,T(-4$ - 2T' - 2T-iT) + Td^{2B' + L" - T') 

+ 2nTd,{2B + L' -T)Y (A.4) 

Se = (AA)"' Qs'S-' - l^S''^ {dj {2B + L'-T) d^T 

+ d^d'^L [2dkdjL + AdkdjE + AEkj + {2nT - 4«')4j] 
+ Td,dj{L' + 2nL) 

+ T {2E[^ + 29,9^^;' + AHE.j + md,djE) 

+ d^Ldk [d.djL + 2^;,^ + 2d^djE) } , (A.5) 

\did^L [2dkdjL + AdkdjE + AEkj + (27ir - 4^')(5fej] 

+ Td^dj (L' + 2HL) + dj {2B + L' ^T) diT 
+ T {2E[j + 2^^^JE' + AUEi^ + AHd^djE) 

+ d'^Ldk {d.djL + 2£;y + 2d^djE) }, (A.6) 

5^ ^ T(p"T + pT' + 2V) + d'Ld^{25p + p'T), (A.7) 

Sp = T{P"T + P'T' + 25P') + d'Ldi{2SP + P'T), (A.8) 

nTd\L' - 2v) + Ta*(2i;' - L") 

+ L^djd\2v - L') + [HT + T' + 2$) 

+ 9^(L'-2t>)9ja*Ll. (A.9) 



Appendix B. Transformation rules of the tetrad fields 



Rab and Sab are defined in Eq. ((52)) . The perturbation variables of the metric are 
defined in Eq. H]). 

Appendix B.0.1. First order 

We can read directly from these expressions the transformation rules for the tetrad 

= r (4^)6^) = - r(a>(i))et - e-;^,a'''r(i?(i)) (b.2) 



=T(vI/(i))eX -eX5'^'=a,.r(i?(i)). 
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Appendix B.0.2. Second order 

rI^J ^ = $(2) - 3$2 ^ g^BQ'B (B.3) 





oo,x 

(2) 

■Oai,X 



m 

aiO,X 

m 

aiak,X 



c(2) _ 
•~'aiO,X - 

c(2) 

aiak ,X 



-s. 



(2) 

oo,x 



^(2) 



- {da.da^Ei^^ + Bill) +i^'5a.a, 
+ 3 d'^^ E + i?,",' ) {da, da, E + Ea,a, ) 

$(2) _ $2 _^ a-sa^B 



In these formulas, we have omitted the first order superscript as there is no possible 
confusion. In the following, we will also omit the first order superscript. The 
transformations rules for the tetrads can be read, as we did for the first order case: 



T 



T 



- [r($(2)) _ 3T($)2 + da{B)d'T{B)\ (B.4) 

+ { - a«'T(S(2)) + [2T($) - 4T(*)] d''^T{B) 

+ 45«^T(B) [d'^*da,T{E)+E'^i^ 
T(*(2)) + 37-(^)2l 

I- { - d'"'da,T{E(^^'>) + 3 [da^d^^TiE) + E^] d^'da^TiE) + E^) 
-6T(*) [5'"=a„,T(i?)+i?«4] }e; 



5/* 



Appendix C. Transformation of (5^^^/ 



(C.l) 



+ 



(1)6 



+ 2 



(1)6 
,X 



These individual terms are explicitly given by 
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r 



(C.2) 

2 



- (t(2)' + HT(2) + 54,(r, L)) + 4$(T' + HT) + 3(T' + HT) 
- 2diBd\-T + L') - di{-T + L')d\-T + L')] 71°^, 



r 



o(l)Oo(l)0 „(1)0^ 



TT TT 



9(7rO)^ 



(C.3) 
(C.4) 



_9 



(T' + WT) ($ + T' + WT) 



(C.5) 



2 ^(7r°)2 {n'diT) ($ + T' + HT) - 2T— ^(c& + T' + UT) 

d{-K^) drjoir 

2^n°(B' - d'T + d'L')diT 

2^7r° [n^d'dj{E + L) + n^E'^ + n^(-* + UT)] diT, 



d 



(C.6) 



r)2 f f)f 



+ 2 



2(7r j + TT-^TT 



<E> (T' + WT + n'^iT) , 



07? aTr'J 



2 , df 



drf 



tT 



dr] 



(TT' + diTd'L) + 



9y~ 



TT°2T{T' + nT + n'diT) 



(C.7) 
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+ ^—{Tiy [2n'^^T{nT + T') + {n'^^T){n^^JT) + H^T'^ + 2nTT' + {T'f] 



5/ 



TT" + H'T^ + mTT' + Tn'd.T' + T'n'^^T + H^T^ + 2nTn'd^T 
+ djT'd^L + djTd^L' + HdjTd^L + d^Ln'didjT + d^Tn'd.djL + [T'f 



2C 



Til 



4^V(a;^ap'')] = (C.8) 

ox^ \ orj ott'J ott* / 

In the above formulas, we have omitted to write the fact that the derivatives with 
respect to rj or are taken at fixed 7r°. 

Appendix D. Integral relations necessary to derive the fluid limit 

The integrations on angular directions can be handled with the general formulas (see 
Ref. [13]) 

d^n 

n'K..n"'—— = Q if fc = 2p + 1 (D.l) 
47r 

/n^i^.n'*— = f(5('i^^..<5*<'=-i)*'=)') if k = 2p. (D.2) 

J 47r K + 1 V / 

By successive integration by parts, we also obtain the following useful results 



/(a;^,7r°)(^°)3d^°d2r! = p(x^). 



Appendix E. The fluid Hmit for radiation 

As explained in section [^TTl second order quantities involve either purely second order 
perturbation variables or terms quadratic in first order perturbation variables. As 
long as the order of the quantity is known we can omit the order superscript in order 
to simplify notations. 

Appendix E.l. Geometric quantities 

In the Newtonian gauge, ignoring vector perturbations for simplicity, the non- 
vanishing Christoffel symbols associated with the metric ([1]) are for the background 

(")r°o = (°)ro, = m,k. ^^^n, = m] . (e.i) 
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At first order, we get 

(i)r[;o = $', (Dro^. = 9,a>, ^^)vi, = d^^, (e.2) 

^^^T% - '2-'HEjk + E'^k - (2H$ + + 5jk, (E.3) 

^^^Tl^ = E" , - ^' 5], (E.4) 

(1) r;., = 29(,[i?;.) - - d\E,k - (e.5) 

where = (^y + Aji)l2. At second order, we obtain 

(^^rfJo = $'-4$$', (2)ro^. ^aj$-4$aj$, (e.6) 

(2) rj,o = - AE'^dj<^ + 4*9*$, (E.7) 

+ 2nEjk - s<PnE,k + s^fe - 4$s;.fc, (e.s) 

-4£;*'=£;[.j- +4*£:;', (e.9) 

(2)r}, = 2%[i?^) - - d^iE.k - ^6,k) (E.IO) 



+ 4(E'^ - *(5' 



Appendix E.2. The radiation fluid equations 

The conservation equation V^T'*" for the stress energy tensor ([7|) with a radiation 
equation of state P = p/3 and where we assume tt^" — 0, are the conservation 
equation 

((5(iV)' + 4H(5(iV+ ^p(a«^^^ -S^-'i'') ==0, (E.ll) 

and the Euler equation 

+ $(1) + ^ ^ 0, (E.12) 
4p 

+ $(2) + ^ ^ 

4p 

where the source terms in the second order equations, which are quadratic in first 
order perturbation variables, are given by 

Sc -|(5p*' + 6p«'*'- ($ + (5)pA?; 
3 V 

+ div [~d'Sp - 2pd'v' - 2p9*$ + 3^9'*] }, (E.13) 
diSe^ ~2(^^div^ +10^'diV + 'i-^d,v' -2dj{d^vd,v) 

+ 2<Pdiv' -2—d,'t> + A<i>d,^. (E.14) 
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